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Abstract 


The breeding method has been implemented in the NASA Seasonal-to-Interannual 
Prediction Project (NSIPP) Coupled General Circulation Model (CGCM) with the goal of 
improving operational seasonal to interannual climate predictions through ensemble 
forecasting and data assimilation. The coupled instability as captured by the breeding 
method is the first attempt to isolate the evolving ENSO instability and its corresponding 
global atmospheric response in a fully coupled ocean-atmosphere GCM. Our results 
show that the growth rate of the coupled bred vectors (BV) peaks at about 3 months 
before a background ENSO event. The dominant growing BV modes are reminiscent of 
the background ENSO anomalies and show a strong tropical response with 
wind/SST/thermocline interrelated in a manner similar to the background ENSO mode. 
They exhibit larger amplitudes in the eastern tropical Pacific, reflecting the natural 
dynamical sensitivity associated with the presence of the shallow thermocline. 

Moreover, the extratropical perturbations associated with these coupled BV modes reveal 
the variations related to the atmospheric teleconnection patterns associated with 
background ENSO variability, e.g. over the North Pacific and North America. A similar 
experiment was carried out with the NCEP/CFS03 CGCM. Comparisons between bred 
vectors from the NSIPP CGCM and NCEP/CFS03 CGCM demonstrate the robustness of 
the results. 

Our results strongly suggest that the breeding method can serve as a natural filter to 
identify the slowly varying, coupled instabilities in a coupled GCM, which can be used to 
construct ensemble perturbations for ensemble forecasts and to estimate the coupled 
background error covariance for coupled data assimilation. 


2 



1. Introduction 


Of all the seasonal-interanmial climate variabilities, the El Nino-Southern Oscillation 
(ENSO) phenomenon plays one of the most important roles in dominating interannual sea 
surface temperature (SST) variability in the tropical Pacific. Feedbacks through strong 
atmosphere-ocean coupling in tropics characterize the co- variability of wind, SST and 
thermocline (or warm water volume) of ENSO, which induces not only a global impact in 
climate anomalies but also modifies the frequency of extreme weather events like floods 
or hurricanes. The delayed oscillator mechanism related to tropical waves propagating 
downwelling/upwelling information in the upper ocean (Schopf and Suarez 1987; Suarez 
and Schopf 1988; Battisti 1988) can explain many features of ENSO. Jin (1997) further 
emphasized the importance of the variations of warm water volume in the upper ocean 
with a warm water recharge/ discharge mech anis m It has been shown that the coupled 
dynamic/thermodynamic mechanisms, i.e., the thermocline and Ekman feedbacks, 
v responsible for delayed or recharge/discharge oscillators, can explain both the west-east 
asymmetry in the climate mean state and the ENSO variability in the equatorial Pacific 
basin (Cai 1995; Jin 1996; Dijkstra and Neelin 1999; Van der Vaart et al. 2000; and Cai 
2003). Those studies point to the essential role of oceanic memory associated with 
oceanic wave dynamics. Through the SST that serves as the lower boundary of 
atmosphere, the oceanic memory dominates the atmospheric seasonal-interannual 
variability. Therefore, being able to estimate seasonal-interannual variations of SST 
becomes a crucial factor for successfully predicting seasonal-interannual variability like 
ENSO. In order to accurately describe the SST and its uncertainties, it is necessary to 
consider variations of subsurface temperature since they are intimately related to SST 
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through thermocline feedback and Ekman-layer dynamics. 

In the last two decades, ENSO prediction skill with dynamical models has been 
steadily improving due to the establishment of the tropical observation network, 
especially the TAO/TRITON mooring arrays, and a better representation of air-sea 
interaction processes in coupled ocean-atmosphere models. The Zebiak and Cane model, 
an intermediate coupled ocean-atmosphere model (Zebiak and Cane 1987, hereafter ZC), 
has been widely used for ENSO studies. However, the presence of errors in the initial 
conditions limits the forecast performance (Latif et al. 1998). More sophisticated 
initialization methods have been shown to be important in improving ENSO prediction. 
For example, data assimilation schemes and ensemble forecasts provide some 
information on the uncertainties embedded in the initial conditions (Chen et al. 1995, 
Rosati et al. 1997 and Ballabrera et al. 2000). Recently, Chen et al. (2004) demonstrated 
that the retrospective forecast from the latest version of this model, due to an 
improvement of the data assimilation to the ZC model, can have a good skill up to two- 
year. They also argued that the evolution of El Nino is controlled to a large degree by 
self-sustaining internal dynamics, suggesting that model-based predictions of El Nino 
depend more on the oceanic initial conditions rather than on unpredictable atmospheric 
noise. 

Several operational forecast centers around the world now use a coupled general 
circulation model (CGCM) to forecast ENSO events with an average lead-time around 
six months. Asummaiy of ENSO forecasts from different operational centers can be 
found under http://iri.columbia.edu/climate/ENSO/currentinfo/ . Current ensemble 
forecasts are based on one of two approaches for the initialization process: the “two tier” 
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and the “one tier” configuration (i.e., double or single stage configuration). In the widely 
used two-tier system, a single forecast of SST anomalies is used to generate an ensemble 
of atmospheric forecasts (Bengtsson et al. 1 993). Since it neglects the coupled nature of 
the initial perturbations, this approach does not project the ENSO mode on those 
perturbations. This hardly seems optimal to for seasonal and interannual prediction. The 
one tier or single stage configuration of CGCM introduced by European Centre for 
Medium-Range Weather Forecasts (ECMWF) (Stockdale et al. 1998) generates all the 
ensemble forecasting members via a coupled ocean-atmospheric model in order to have 
the perturbation growing under a coupled configuration. The NSIPP forecast ensemble, 
also a one-tier system, includes perturbations in both the atmosphere and ocean, but the 
initial perturbations are generated independently. Although coupled instabilities will 
eventually develop in these one-tier systems, they still handicapped by not including 
coupled uncertainties in the initial perturbations for ensemble predictions. Therefore, 
there is a need for ensemble ENSO prediction systems to include coupled initial 
perturbations and feedbacks. 

The breeding method (Toth and Kalnay 1993, 1997) and singular vectors (Errico and 
Vukicevic 1992; Buizza and Palmer 1995; and Palmer et al. 1998) are two of the main 
methods used operationally for generating effective ensemble members in ensemble 
forecasts with an atmospheric GCM model. Studies of the growth of error/growing 
modes related with coupled ocean/atmosphere instabilities have generally focused on 
obtaining singular vector (SV) from a linear/adjoint propagator of intermediate models 
(e.g., the ZC model). Chen et al. (1997) showed that the error growth depends on the 
phases of the ENSO and the seasonal cycles, even though initial and final SVs are 
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insensitive to both. It should be pointed out that their analysis was done using the linear 
propagator with SST perturbations only, neglecting the coupled nature of the 
perturbations of other variables. Xue et al. (1997 a,b) used a reduced EOF space 
spanning the ZC solutions in order to be able to compute the SV for all model variables 
in a full tangential linear ZC model. They showed that the growth rate and singular 
vectors are similar to Chen et al. (1997) when choosing an SST norm. However, with an 
energy norm, the wind and thermocline fields become more important than the SST. Fan 
et al. (2000) found that SST plays a very important role in seasonal (3-6 months) 
predictability, and pointed out that using the analysis error covariance as a norm in 
calculating SV yields quite different results than the frequently used energy norm. The 
strong dependence of singular vectors on the choice of norm (i.e., the definition of the 
“size” whose growth is to be maximized) and the choice of optimization time, become 
limitations when applying this method to a complex system like CGCMs because of 
existence of various types of instabilities over a whole spectrum of scales. It is difficult 
to cleanly separate these modes and to keep the coupled instability as the dominant 
growing mode, since the adjoint model, being linear, can be swamped by the presence of 
much faster growing atmospheric and oceanic instabilities (Pena and Kalnay, 2004). 

Toth and Kalnay (1996) suggested applying the breeding method in a coupled 
ocean-atmosphere system to isolate ENSO coupled instabilities. The bred perturbations 
are a superposition of the leading Lyapunov vectors (dominant instabilities) in the 
dynamical system and the advantages of this method are that it is simple, efficient and 
independent of the choice of norm. Cai et al. (2003) first tested the breeding method in a 
coupled system using the ZC model. They found that bred vectors are capable of 
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describing the characteristics of coupled instabilities associated with the ENSO 
development. They found that the growth rate is weakest at the peak time of the ENSO 
states (both positive and negative) and strongest between the events. Unlike the singular 
vectors, die coupled bred vectors are insensitive to the choice of norm and very sensitive 
to the background ENSO phase and the time of the year. Their results also suggest that 
the presence of the “spring-barrier” in the ENSO prediction skill is related to the coupled 
instability. Filtering out the coupled bred mode from initial perturbations greatly 
increases the time for doubling forecast error and reduces the “spring-barrier”. They also 
found that using a pair of plus/minus bred vectors as ensemble perturbations led to a 
significant improvement in the ensemble mean forecast. These results illustrate the 
potential impact of coupled bred vectors in both ensemble prediction and data 
assimilation for ENSO predictions. Pena and Kalnay (2004) tested the breeding method 
' for coupled Lorenz models with distinct time scales to mimic the interaction between a 
slow “coupled tropical ocean-atmosphere system” and an “extratropical atmosphere”. 
They found that breeding is able to isolate the slow modes of the coupled system when 
rescaling intervals and amplitudes are chosen from physically appropriate scales and the 
rescaling factor is obtained from the slow component of the system. In contrast, 
Lyapunov and Singular vectors are unable to isolate the slow modes, because they are 
based on linear models, and are therefore dominated by the fast modes. 

The results with simpler models encourage us to implement the breeding method in 
a more complicated and complete model, like a CGCM, that includes many types of 
instabilities, without sacrificing resolution or simplifying model physics. As a first step, 
we examine in this paper the growing coupled instabilities obtained using breeding in a 
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fully coupled general circulation model. Our objective is to identify the characteristics of 
bred vectors associated with the ENSO mode derived from the CGCM and to investigate 
whether an ENSO-related coupled instability can be isolated from weather noise using 
the breeding method. Specifically, in this paper, we would like to address the following 
questions: (1) Can breeding be used to identify the coupled, slowly growing.ENSO 
instability and isolate it from other short term, small scale instabilities such as weather 
noise? (2) How does the coupled instability evolve with the background ENSO 
variability? and (3) Are the main characteristics of coupled bred vectors reproducible 
with two different coupled ocean-atmosphere GCMs? 

The structure of the paper is as follows. In section 2, we give a brief description of 
the NSIPP coupled model, which has been used to generate coupled bred vectors. A brief 
discussion about the ENSO variability in the NSIPP model is also included in section 2. 
Section 3 describes how the breeding method is applied in a coupled GCM. Section 4 is 
devoted to describe the main characteristics of coupled bred vectors derived in the NSIPP 
model. A comparison of the results obtained from NSIPP and from NCEP/CFS03 is also 
presented in section 4. A brief summary and discussion of the next phase of our research 
are included in section 5. 

2. The NSIPP coupled global circulation model and its ENSO variability 

In this study, we test the breeding method on the NSIPP coupled ocean-atmosphere 
general circulation model, in a perfect model scenario. The NSIPP coupled model is a 
fully coupled global ocean-atmosphere-land system developed at NASA Goddard Space 
Flight Center (GSFC) (Miller et al. 2004; Vintzileos et al. 2003). It is comprised of the 
NSIPP-atmospheric model (AGCM, Bacmeister and Suarez 2003; Bacmeister et al. 
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2000), the Poseidon ocean model (OGCM, Schopf and Loughe 1995), and the Mosaic 
land surface model (LSM, Koster and Suarez 1992). The NSIPP AGCM uses a finite- 
difference C grid in the horizontal and a generalized sigma coordinate in the vertical. An 
important feature of this AGCM is its inclusion of an empirical cloud diagnostic model 
and a relaxed Arakawa-Schubert cumulus/boundary-layer parameterization. The OGCM 
is a Poseidon quasi-isopycnal model designed with generalized horizontal and vertical 
coordinates. The current NSIPP version uses a reduced-gravity formulation and an 
embedded surface mixed layer following Sterl and Kaltenberg (1994). The isopycnal 
region is treated in a quasi-isopycnal fashion, in which layers do not vanish at outcrops 
and remain a thin minimum thickness at all grid points. The treatment of horizontal 
mixing within the model is implemented with high order Shapiro filtering. The NSIPP 
CGCM employs the Goddard Earth Modeling System (GEMS) to couple the atmosphere, 
ocean and land models. The ocean and atmosphere exchange information once a day. 

- Operational forecasts and hindcasts indicate that forecasts for the Nino3 index 1 remain 
skillful up to 9 months lead-time, depending on the starting month. Forecast information 
can be found under http://nsipp.gsfc.nasa.gov/ . 

A 50-year simulation run has been made with a research version of the NSIPP 

CGCM identical to the operational model except for a slightly coarser resolution 

(AGCM: 3.75° in longitude by 3° in latitude and 34 sigma layers; OGCM: 1.25° in 

longitude by 1/2° in latitude and 27 layers). This integration is referred to as the “control” 

or “background” run in this paper. The Nino3 index obtained from this control run (not 

shown) exhibits a realistic ENSO-like variability, although its biennial component is 

1 The Nino3 index is defined as the spatial average of SST anomalies (control) in the Nino3 domain (150 E- 
90°W; 5°S-5°N). The BV Nino3 index is defined in the same way except that the B V SST perturbations are 
used instead of SST anomalies. 
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stronger than observed. In agreement with observations (Boulanger and Menkes 1995, 
1999), the control run simulation shows warm anomalies at the thermocline level prior to 
the SST anomalies. This indicates that air-sea interaction processes and coupled 
instabilities can be reasonably represented by the CGCM. The spatial patterns of the 
coupled ENSO variability of the control run are presented in the left panels of Figs. 4-6 
and will be discussed together with the presentation of the bred vector spatial structure in 
Section 4. 

3. Breeding method with the coupled NSIPP CGCM 

We used 10 years of the control run, from model year 2020 to 2029, to perform 
breeding experiments (discussed in the next section) under a “perfect model scenario”, 
i.e., assuming that the control run is the “truth” (or analysis in a forecast system). The 
breeding method (Toth and Kalnay 1993, 1996, and 1997) originally developed to 
“breed” fast growing modes for short-to-medium-range atmospheric ensemble 
forecasting was implemented in the NSIPP CGCM. The breeding cycle includes the 
following steps (Fig. l): (1) add a random perturbation on the initial analysis fields (both 
oceanic and atmospheric restart files); (2) make a one month forecast (coupled breeding 
run) starting from the perturbed fields; (3) take the difference between the breeding and 
the control runs; (4) rescale the difference field (referred to as the “bred vector”) to the 
same size as the initial perturbation; (5) add the rescaled bred vector field to the next 
analysis field; and (6) repeat steps (2) to (5), adding the coupled bred perturbations onto 
analysis fields and evolving them with background flow, which forms the “breeding 
cycle”, throughout the 10 years. Thus, bred vectors are perturbations aligning along 
favorable growing directions, periodically resized and added to the next analysis field in 
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each of the breeding cycles. Bred vectors are essentially a finite-time, non-linear 
extension of local Lyapunov vectors, representing the preferred instability directions of 
the evolving basic flow. Under the “perfect model scenario”, the analysis fields are 
considered to be the truth so that both analysis and forecast fields refer to the same 
control run. 

When the system contains many different types of instabilities, like the coupled 
ocean-atmosphere model, the initial (or rescaled) perturbation amplitude and the 
rescaling time interval are the two parameters that can be used to select the type of the 
instability that shows up in the bred vectors (Toth and Kalnay 1993, 1997; Pena and 
Kalnay 2004). The BVs consist of perturbations related to instabilities whose saturation 
amplitude is well above the range of amplitudes allowed, and whose growth rate is 
largest, since perturbations related to slower-growing instabilities are damped more 
strongly by the rescaling cycle. For example, convective instabilities will quickly 
saturate at very small amplitude when breeding in a global atmosphere. A similar 
situation appears in the instability of a coupled atmosphere-ocean system, but in this case 
the coupled atmospheric signal that we are seeking has smaller amplitudes than the 
atmospheric “weather noise”. Therefore, choosing an atmospheric measure of the 
amplitude for the rescaling would not be an effective way of retaining the slower growing 
coupled perturbations (Pena and Kalnay, 2004). For this reason, the variables we choose 
to rescale the perturbations should be primarily based on oceanic quantities whose 
variability is dominant at the seasonal to interannual time scales. Moreover, the rescaling 
time interval in a breeding cycle plays a crucial role for capturing coupled instabilities 
corresponding to seasonal to interannual time scales and for isolating them from weather 
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noise. Choosing a rescaling time interval of one month allows saturation of the fast 
growing weather noise. 

We implemented the breeding method in the NSIPP global CGCM by choosing the 
period of rescaling as one month and a Nino3 SST root-mean-square norm (RMS), as in 
Cai et al. (2003) in the Zebiak and Cane model. We chose the RMS SST amplitude 
within the tropical Pacific domain (15°S-15°N and 120°E-90°W) of 0.085°C, or about 
10% of the background SST variability. In the perfect model configuration, the control 
run is used as both the analysis and the background (first guess). The breeding cycle is 
simplified by adding the perturbations to the control run, rather than to an analysis, and 
the bred vectors are the rescaled difference between the perturbed run and the control run. 
Two independent breeding runs were made starting from two independent random 
perturbations, which were created by taking the difference between two model states at 
randomly chosen months. Each run contains 123 months, starting from model year 2019 
September to year 2029 December. The starting month is chosen so that the first major 
warming event takes place about 2 years into the breeding run to ensure the coupled bred 
modes are closely associated with ENSO. The first three months of the breeding runs are 
treated as a transient period allowing the bred vectors to grow from random perturbations. 
The analysis presented below is derived from the remaining 120 breeding cycles (120 
months). We found that the two independent 10-year breeding cycles yielded very 
similar bred vectors (not shown), so that we combine the two bred vector perturbations as 
a single time series of 20 years to reduce sampling errors. Hereafter, we will refer to the 
combined bred vectors as BV perturbations. 

4. Bred vectors in the NSIPP coupled model 
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4.1. Growth rate of coupled bred vectors 


Bred vectors represent, by construction, the instabilities that grow on the background 
flow. An example of the instabilities captured by the coupled bred vectors is shown in 
Fig. 2, a snapshot of the bred vector SST field (contours) together with the corresponding 
background SST field (shading) on July 1 of the model year 2024. It shows that the bred 
vector field has large amplitude along the sharp temperature gradient in the equatorial 
cold tongue, coinciding with the background waves along the edge of the cold tongue. 

The bred perturbation seems to suggest that the instability tends to make the waves break. 
Clearly, the formation of tropical instability waves is captured by the bred vectors. 
Although the bred vectors may include both coupled and uncoupled ocean instabilities 
such as those shown in Fig. 2, we illustrate in the remaining part of the paper that a 
significant portion of the instabilities are related to coupled ENSO dynamics. 

The growth rate of the coupled bred vectors is calculated based on the chosen 
rescaling norm of the perturbation field within the tropical Pacific region: 


G(t) = 


J - 

V NG 

J- 

Y NG 


( 1 ) 


where NG is the total number of model grid points in the Nino3 region and t is the model 
time in months. In other words, we measure the growth rate of bred vectors by their 
amplification factor within a month. In the 10-year experiment, the typical value of the 
growth rate in one month of the NSIPP coupled model varies around 3 to 5 (not shown), 
which is much larger than the coupled instability found by Cai et al. (2003) for the ENSO 
mode in the ZC model, about 1 to 3 per month. This is to be expected since the growth 
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rate in the coupled GCM includes both coupled and uncoupled instabilities of any kind, 
such as the tropical wave instability shown in Fig. 2. We interpret the growth rate of 3-5 
per month as consisting of a noisy background growth rate of (mostly uncoupled) 
instabilities of about 3 per month (a component essentially absent in the ZC model), plus 
a coupled growth of about 1-2 per month, associated with the ENSO signal that we are 
seeking. 

In order to test whether there is a component of the perturbation growth (above the 
background noisy growth rate of about 3 per month) evolving upon the coupled ENSO 
background state (rather than growing randomly), we calculate the lag/lead correlation 
between the growth rate and the absolute value of the background Nino3 index. We use 
the absolute value of the Nino3 index in order to account for the large amplitude of both 
positive and negative SST anomalies. It is evident in Fig. 3 that the growth rate of 
coupled bred vectors tends to be largest about 3-4 months prior to the time at which the 
background ENSO amplitude reaches its maximum stage (positive or negative). There is 
also a relatively weaker growth at about 4 months after the maximum stage of the 
background ENSO. These results are qualitatively in good agreement with the results 
obtained with the ZC model in Cai et al. (2003). The results shown here suggest that the 
breeding method can serve as a natural filter to identify the slow, coupled instability 
related to ENSO variability by selecting the proper rescaling parameter in the breeding 
cycle. Therefore, we can expect that projecting the initial perturbations of ensemble 
members onto ENSO-related growing errors should improve predictions in an ensemble 
forecast system, as in Cai et al (2003). 

4.2. The structure of the coupled BV mode 
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The spatial patterns of the coupled BY can be identified by constructing regression 
maps for both oceanic and atmospheric variables against the BV Nino3 index, defined as 
the spatial average of BV SST in the Nino3 domain. This regression filters out 
perturbations unrelated with ENSO, like the tropical equatorial instability shown in Fig.2. 
We will compare those maps with the background regression maps constructed with the 
same regression method but using the background Nino3 index in order to determine 
whether the coupled BV modes can capture the growing features associated with 
background ENSO variability. 

The oceanic global regression maps for the background fields show the typical 
tropical variability corresponding to the ENSO mature stage (Fig. 4(a)-(c)). These 
' patterns include a large warming extended from the east to central equatorial Pacific, a 
deepening thermocline in the eastern equatorial Pacific, and an accompanying shallowing 
feature off tire equator in the western basin, and a basin-wide eastward current anomaly. 

- The regression maps for the BV fields are shown in Fig. 4(d)-(f). The coupled BV mode 
exhibits a strong signal in the equatorial Pacific and fairly weak variability away from the 
tropics. The patterns of the coupled B V mode are reminiscent of those in the background 
state except that the BV mode is more confined to the east and to the equator. This 
feature is physically meaningful, since it reflects a larger sensitivity to perturbations of 
the background flow in the shallower thermocline in the east along the equator. It is also 
consistent with the delayed oscillator theory, which considers that the perturbations grow 
primarily over the eastern equatorial basin. It is known that in the mean the easterly wind 
stress is balanced by the zonal pressure gradient, which piles up warm water in the 
western basin, resulting in a thermocline that slopes down toward the west. The shoaling 
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thermocline in the east implies that the thermodynamic feedback between SST and near 
surface ocean variables is much stronger in the east than in the west. As a result, oceanic 
perturbations in the eastern basin will be easily amplified through positive feedbacks 
from air-sea interaction. 

The atmospheric components of the ENSO mode derived from the control run and 
the BV field are displayed in Fig. 5. We first examine the tropical Pacific domain to see 
the direct impact related to boundary heating anomalies. Fig. 5 shows that the patterns in 
the BV fields have many features in common with the patterns of the background state, 
namely, westerly wind perturbations located in the central equatorial Pacific and the 
baroclinic structure in the height fields corresponding to the location of BV SST in Fig 
4(d). In addition, the BV outgoing radiation reflects an enhanced convection activity in 
the eastern basin. These features reinforce the conclusion that the leading coupled BV 
mode is related to the coupled instability. 

It is of interest to point out that the coupled BV also reflects the sensitivities in 
extratropical regions associated with background ENSO atmospheric teleconnections. 
Shown in Fig. 6 are the regression maps of surface pressure and geopotential at 200mb in 
Northern Hemisphere for the background state and for the BV field. The teleconnection 
patterns of the background state indicate a low-pressure anomaly over the North Pacific 
and a high-pressure anomaly over North America. This barotropic structure is very 
robust and extends to a high altitude. It is induced by wave-train patterns associated with 
the large scale heating in the tropics. For BV maps, strong responses can also be 
identified in those regions, especially where background regression maps show a strong 
gradient, for example in the mid Pacific at 30°N and east coast of North America. Wave- 
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train patterns can also be found in BV regression maps. In the Southern extratropical 
region (not shown), atmospheric regression maps show a teleconnected pattern associated 
with background ENSO, and related BV dynamical sensitivities. 

In order to show how the coupled BV mode evolves with the background ENSO 
evolution as in Cai et al. (2003) with the ZC model, we construct lead/lag regression 
maps against the time series of the amplitude of the background Nino3 index with a 
lead/lag time up to 6 months (Fig. 7). It is clear that the temporal evolution of the 
coupled BV mode is highly related to the background ENSO evolution. It shows that in 
the eastern basin the coupled BV mode leads the large amplitude of the background 
ENSO events by several months. The signal is clearly coming from the coupled 
dynamics because an increase of the ocean heat content and a warm SST anomaly in the 
eastern basin as well as the presence of westerly wind anomalies all lead the background 
ENSO events by about 2-3 months. This coincides with the timin g of the maximum 
growth of the coupled BV mode which also leads the background ENSO events by 3 
months (Fig. 3). It is seen from Fig. 7 that west of 130°W, BV surface height and zonal 
wind stress exhibits a lag response of the background ENSO of about 3 months. 

Therefore, using the breeding method, we show that the dominant instability in the 
NSIPP GCM model is initiated in the eastern basin and that the signal is a coupled 
instability. 

4.3. Comparison of the NSIPP and the NCEP/CFS03 Bred Vectors 

Similar breeding experiments were carried out with the coupled forecast system 
model (CFS03) developed in the National Centers for Environmental Prediction (NCEP). 
The atmospheric component uses the current version of medium range forecast (MRF) 
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global model with a spectral truncation of 62 waves (T62) in the horizontal (equivalent to 
nearly 200 Km) and 64 vertical levels in sigma coordinate (Kanamitsu 1989; Kanamitsu 
et al. 1991; Caplan et al. 1997; Wu et al. 1997). The ocean component is the GFDL 
Modular Ocean Model V.3 (MOM3) with 40 layers in the vertical (Pacanowski and 
Griffies, 1998). The zonal resolution is 1° and the meridional resolution is 1/3° between 
10°S and KPN, gradually increasing through the tropics until it is fixed at 1° poleward of 
30°S and 30°N. 

Two independent breeding experiments were performed by choosing the last 4 years 
from a 23-year perfect model experiment as the background state. This 4-year period 
covers a warm event which matures at the model year 21,2 years into the breeding run. 
The rescaling factor for perturbations is based on the SST norm in the whole tropical belt 
(10S-10N) and the perturbation size was chosen as 0.1°C. As in the breeding 

experiments performed with the NSIPP CGCM, we chose one-month as the rescaling 
period. Like the NSIPP coupled experiments, the two BV runs for the NCEP system 
were very similar despite having been started with different random perturbations so that 
their results are processed as a single 8-year time series. Comparisons between the 
results from the NSIPP and the NCEP/CFS03 coupled system are made for the purpose of 
demonstrating the robustness of the bred vectors in coupled GCMs 2 . 

Fig 8 (a)-(f) are background oceanic regression maps of two coupled GCMs. The 

oceanic components from the two GCMs successfully produce fundamental features of 

ENSO. Their differences also reflect differences of numerical schemes in the model 

dynamics or different choices of physical parameterizations. The meridional structure of 

2 Unfortunately, the experiments performed at NCEP were erased, so that we have only a limited number of 
diagnostic comparisons available. 
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warming and thickening surface height (Figs. 8(a)-(b)) in the NCEP/CFS03 GCM in the 
Eastern Pacific is wider than that of the NSIPP GCM (Figs. 8(d) and (e)). In addition, the 
regressed surface height of NCEP/CFS03 shows the southern branch of the shoaling 
patterns off the equator extends more southward instead of being meridionally limited as 
in the NSIPP case. This can also be seen in the SST and zonal current patterns. Despite 
those distinctions, the bred vectors from the two coupled systems have significant 
similarities linked with the background ENSO. To compare BV structures, we show the 
EOF modes of oceanic variables. Fig. 9(a)-(c) are the first EOF mode of the BV SST and 
first two modes for the BV Ihermocline from the NSIPP CGCM and Fig. 9(d)-(f) are the 
same modes using BV from NCEP/CFS03. Despite the fact that these are two different 
CGCMs, with very different background evolution, there is a strong resemblance between 
the BV EOF modes. Both die leading modes (EOF1) in NSIPP and NCEP/CFS03 bred 
vectors based on SST show an ENSO-associated warm feature in the tropical Eastern 
Pacific, farther east than in their respective background (not shown). Reflecting the 
different mean structures and background ENSO variabilities from a different coupled 
system, the NCEP/CFS03 BV SST EOF1 extends over a larger spatial scale, covering the 
whole Nino3 domain while the BV SST EOF 1 from NSIPP model is confined to east of 
130°W and is meridionally limited. These two EOF1 modes respectively explain 11% 
and 14% variance from the total growing SST perturbations. This suggests that the 
coupled growing perturbations associated with ENSO variability represent at least 10% 
of die total growing perturbations due to a variety of instabilities that appear in a coupled 
GCM. The fact that the leading EOF modes from BV fields in both coupled systems show 
an ENSO-like structure confirms our conjecture that the breeding method is capable of 
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capturing the coupled instability even in the presence of other types of instabilities in the 
fully coupled GCM model. Moreover, this mode is robust and dominant. Even if we 
enlarge the domain for the EOF analysis to a global domain, the same mode appears, 
showing fairly weak oceanic amplitude in the extra-tropics. This indicates that the 
breeding method can help to identify the largest growing error projecting on the ENSO 
variability in a global coupled model and can be applied to a model with full, complete 
physics such as a GCM, with the full resolution of model output. 

Similar natural sensitivities in the eastern Pacific can also be found in the B V 
thermociine fields of both systems (Figs. 9 (b), (c), (e) and (f)). Both leading EOF modes 
of the BV thermociine have a deepening feature along the equator and shoaling features 
off the equator, except that the EOF1 from NCEP/CFS03 extends almost across the whole 
basin. In addition, both EOF2 modes have a dipole pattern along the equator and 
establish a wave couplet off the equator in the western basin. All the information above 
reveals that oceanic perturbations will develop as Kelvin/Rossby wave packages, 
propagating the upwelling/downwelling signals in the tropical region. It should be noted 
that the EOFl of the BV thermociine is close to the EOFl of BV SST and there is a high 
correlation between their corresponding leading principal components (not shown). 

These two modes represent the dominant growing coupled instability, which has been 
also been obtained by the BV oceanic regression maps. The robustness of the results 
from two different coupled models supports our hypothesis that bred vectors are 
associated with the background ENSO variability. The differences between the BVs 
indicate that bred vectors are sensitive to model behavior. For example, different vertical 
mixing schemes adopted in ocean models will have an impact on thermociine variations 
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particularly in the shallow mixed layer region. 

Based on the regression maps against with BV Nino3 index (not shown) in the 
tropical region, we can also deduce that the coupling strengths are different in the two 
coupled GCMs. In both ocean components, a 1 -meter variation (deepening) in BV 
thermocline corresponds to 0. 1 °C warming in the eastern Pacific. The corresponding BV 
zonal wind stress shows a perturbation of 1.5 Nm~ 2 from the NSIPP AGCM, and it 
prevails in the central basin. In contrast, the corresponding stress perturbation is only 0.5 
Nm’ 2 in NCEP/CFS03 case. In addition, the regressed BV surface pressure and 
geopotential height for the NCEP/CFS03 model are less organized in the tropics than for 
the NSIPP coupled model. This seems to suggest that perturbations are more strongly 
coupled in the NSIPP CGCM in the tropical domain. 

There are also similarities between the two systems in the extratropical ENSO- 
associated teleconnection patterns. Fig. 10(a),(b) are the regression maps of BV surface 
pressure from the two coupled models. Similar responses can be identified from the 
eastern basin of the North Pacific to the North Atlantic despite the different responses for 
other locations. This type of information should be particularly valuable for ensemble 
forecasting since the robust feature teleconnected to ENSO variability will remain and 
influence other locations. 

5. Summary 

In this study, we demonstrated for the first time the feasibility of applying the 
breeding method to a global coupled ocean-atmosphere general circulation model to 
obtain ENSO coupled instabilities. The results from the BV fields derived from the 
NSIPP coupled model show that the breeding method is capable of obtaining die fast 
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growing coupled modes associated with ENSO variability. Potential predictive impact 
has been shown in the BV growth rate, by the fact that the maximum growth leads the 
warm/cold event by about three months. Therefore, the growth rate can be viewed as a 
precursor to background ENSO variability. In particular, the amplitude of the BV in the 
eastern tropical Pacific increases before the development of the background El Nino 
events. Therefore, the bred vectors contain ENSO-related perturbations that can be used 
as the initial perturbations for ensemble forecasting for improving ENSO forecast skill. 
The three-month lead time in our results (compared with 6-12 months for the ZC model) 
may seem to be too short for improving current ensemble forecast system. However, we 
conjecture that this relatively short t im e-lead is partially due to the fact that the ENSO 
cycle in this long simulation is more biennial than 3-7 years. This biennial time scale 
shortens the “growing season” of bred vectors because of the relatively fast pace of the 
evolving background state. We would expect to obtain a longer lead time from a BV 
growth rate derived from the NSIPP operational system since it will reflect the real 
oceanic memory from the observations during initialization, although the shorter time-lag 
may be also associated with the presence of other weather and oceanic instabilities. 

One of the main tasks in our study has been to extract the physical growing mode 
from the total growing variability, separating its signal from noise associated with 
weather and other instabilities. We used regression BV maps to illustrate that the ENSO- 
associated features can be captured by the breeding method. For the oceanic fields, the 
variability mainly comes from the tropical Pacific and has many features in common with 
the background ENSO variability. The mean structure makes the eastern equatorial 
Pacific, with a shallow thermocline, the region of strongest dynamical sensitivity. 
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Corresponding to the BV oceanic fields, the atmospheric BV response is also reminiscent 
of the background ENSO features (e.g. a surface pressure anomaly field with a sea-saw 
pattern from west to east) in the tropical Pacific, indicating that the dominant mode of BV 
is indeed related to ENSO-associated coupled instabilities. In addition to the coupled 
characteristics shown in the tropical domain, the extratropical circulation anomalies 
associated with these coupled BV display a wave-train like teleconnection pattern over 
the North Pacific and North America. The patterns are known to be strongly 
teleconnected to ENSO variability. 

The robustness of the coupled BV modes has been demonstrated by comparing bred 
vectors derived from the NSIPP CGCM and NCEP/CFS03 CGCMs. In an EOF analysis 
applied to BV fields, the leading EOF modes for both systems show similar ENSO- 
related features (warming/deepening thermocline) in the tropical eastern Pacific, 
indicating that the leading fast growing mode is due to the ENSO coupled instability. In 
addition, a strong resemblance between these two independent experiments can be found 
in many fields, even in those atmospheric teleconnected regions associated with 
background ENSO development Our results indicate that the breeding method can serve 
as a natural filter to isolate the slowly-varying, coupled instabilities in a coupled GCM. 
Furthermore, global sensitivities associated with the coupled instability initiated from the 
tropical Pacific can be retained even though the rescaling is simply done in the tropical 
Pacific. 

We conjecture that capturing the coupled bred vectors may benefit data assimilation 
and ensemble forecasting and improve ENSO prediction skill, as shown in Cai et al 
(2003) for the ZC model. Therefore, the next stage of our research is to employ these 
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methods in the NSIPP operational system. In particular, we plan to test whether the 
ensemble forecasts with coupled BVs as initial perturbations are more effective than 
current ensembles that use random perturbations introduced in the separate components. 
Corazza et al (2002) found that in a quasi-geostrophic system, a hybrid 3DVar 
background error covariance augmented with the inner product of BVs improved the 
analysis and the forecasts. Similarly, for the coupled system, the BVs should give the 
structure of the “enrors of the month”, so that these structures could be used to augment 
the background error covariance of the current Optimal Interpolation data assimilation 
scheme and improve the assimilation and forecast of the ENSO-related instabilities. 
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Figure Captions 


Figure 1 Schematic diagram showing a continuously evolving breeding cycle upon 
the unperturbed (control) model integration in a perfect model setting. The difference 
between the unperturbed and perturbed forecasts yields the bred vectors. The growth rate 
is computed as the ratio of the final to the initial size. 

Figure 2 A snapshot of SST in Eastern Pacific showing the bred vectors 
perturbation (contours: Cl = 0. 15 °C) evolving with the background flow (shading with 

an interval of 1°C from 21°C to 30°C) on July 1 of the model year 2024. The dotted 
contours of the BV indicate negative values. 

Figure 3 Lead/lag correlations between the BV growth rate and the absolute value 

of the background Nino3 index. 

Figure 4 Oceanic regression maps in the global domain. Left panels are the 
background fields and right panels the BV fields, (a) SST anomaly (°C); (b) Z20 

anomaly (m); (c) surface zonal current anomaly (ms' 1 ); (d) BV SST (°C); (e) BV Z2 0 

(m); and (f) BV surface zonal current (ms 1 ). Background fields are regressed with the 
background Nino3 index and BV fields regressed with the BV Nino3 index. The scales of 
B V fields are arbitrary but the ratio among B V variables (both oceanic and atmospheric 
variables) is retained as in the original BV fields. 
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Figure 5 The same as Figure 4 except for the atmospheric regression maps over the 
equatorial Pacific basin, (a) wind field anomaly at 850mb (ms' 1 ); (b) surface pressure 
anomaly (mb); (c) geopotential anomaly at 200mb (mV 2 ); (d) outgoing long wave 
radiation (Wm' 2 ); (e) BV wind field at 850mb (ms* 1 ); (f) BV surface pressure (mb); 

(g) BV geopotential at 200mb (mV 2 ); and (h) BV outgoing long wave radiation (Wm' 2 ). 

Figure 6 The same as Figure 4 except for atmospheric regression maps over the 
Pacific portion of the Northern Hemisphere, (a) background surface pressure anomaly 
(mb); (b) background geopotential anomaly at 200mb (mV 2 ); (c) BV surface pressure 
(mb); and (d) BV geopotential at 200mb (mV 2 ). 

Figure 7 Lead/Lag regression maps along the equator for BV oceanic fields against 
the absolute value of the background Nino3 index, (a) SST (°C); (b) surface height (m); 
and (c) zonal wind stress (Nm' 2 ). The contours are arbitrary but the ratio among BV 
variables is retained as in the original BV fields. 

Figure 8 Background oceanic regression maps for two coupled GCMs in the 
tropical Pacific domain. Left panels are the NSIPP anomalies of NSIPP and the right 
panels the NCEP/CFS03 anomalies, (a) NSIPP SST (°C); (b) NSIPP surface height (m); 
(c) NSIPP surface zonal current (ms' 1 ); (d) NCEP SST (°C); (e) NCEP surface height (m); 

and (f) NCEP surface zonal current (ms' 1 ). The regression maps of NSIPP (NCEP) fields 
are computed using the NSIPP (NCEP) Niiio3 index, respectively. 
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Figure 9 The leading EOFs of the BV SST and Z20 perturbations derived from the 

NSIPP and NCEP/CFS03 CGCMs. (a) EOF1 of NSIPP BV SST; (b) EOF1 of NSIPP BV 
Z20; (c) EOF2 of NSIPP BV Z20; (d) EOF1 of NCEP BV SST; (e) EOF1 of NCEP BV 
Z20; and (f) EOF2 of NCEP BV Z20. The scale is arbitrary. 

Figure 10 Atmospheric regression maps of BV 500mb geopotential height in in the 
Northern Hemisphere (m). (a) for NSIPP and (b) for NCEP/CFS03. Both fields are 
computed against their own B V Nino3 indices. 
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Figure 1 Schematic diagram showing a continuously evolving breeding cycle 
upon the unperturbed (control) model integration in a perfect model setting. 
The difference between the unperturbed and perturbed forecasts yields the bred 
vectors. The growth rate is computed as the ratio of the final to the initial size. 
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Figure 2 A snapshot of SST in Eastern Pacific showing the bred vectors 
perturbation (contours: Cl = 0.15 °C) evolving with the background flow 
(shadings with an interval of 1°C from 21°C to 30°C) on July 1 of the model 
year 2024. The dotted contours of the BY indicate negative values. 
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Figure 3 Lead/lag correlations between the BV growth rate and the absolute 
value of the background Nino3 index. 
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Figure 4 Oceanic regression maps in the global domain. Left panels are the 
background fields and right panels the BV fields, (a) SST anomaly (°C); (b) Z20 
anomaly (m); (c) surface zonal current anomaly (ms' 1 ), (d) BV SST (°C); (e) BV 
Z20 (m); and (f) BV surface zonal current (ms* 1 ). Background fields are 
regressed with the background Nino3 index and BV fields regressed with the BV 
Nino3 index. The scales of BV fields are arbitrary but the ratio among BV 
variables (both oceanic and atmospheric variables) is retained as in the original 
BV fields. 
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(c) CNT geopotential anomaly at 200mb 



(d) CNT OLR anomaly 



(g) BV geopotential at 200 mb 



(h) BV OLR 



Figure 5 The same as Figure 4 except for the atmospheric regression maps over the 
equatorial Pacific basin, (a) wind field anomaly at 850mb (ms' 1 ); (b) surface 
pressure anomaly (mb); (c) geopotential anomaly at 200mb (m 2 s‘ 2 ); (d) outgoing 
long wave radiation (Wm‘ ); (e) BV wind field at 850mb (ms' 1 ); (f) BV surface 
pressure (mb); (g) BV geopotential at 200mb (mV 2 ); and (h) BV outgoing long 
wave radiation (Wm' 2 ). 
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(a) CNT Surface Pressure (c) BV Surface Pressure 



Figure 6 The same as Figure 4 except for atmospheric regression maps over the 
f Pacific portion of the Northern Hemisphere, (a) backgroimd surface pressure 
anomaly (mb); (b) background geopotential anomaly at 200mb (m 2 s' 2 ); (c) BV 
• surface pressure (mb); and (d) BV geopotential at 200mb (m 2 s' 2 ). 
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(a) BV SST vs. [Background NIN03| 



(b) BV Hs vs. [Background NIN03[ 



(c) BV zonal windstress vs. [Background NIN03| 



Figure 7 Lead/Lag regression maps along the equator for BV oceanic fields against 
the absolute value of the background Nino3 index, (a) SST (°C); (b) surface height 
(m); and (c) zonal wind stress (Nm‘ 2 ). The contours are arbitrary but the ratio 
among BV variables is retained as in the original BV fields. 
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Figure 8 Background oceanic regression maps for two coupled GCMs in the 
tropical Pacific domain. Left panels are the NSEPP anomalies and the right panels 
the NCEP/CFS03 anomalies, (a) NSEPP SST (°C); (b) NSDPP surface height (m); 
(c) NSIPP surface zonal current (ms 1 ); (d) NCEP SST (°C); (e) NCEP surface 
height (m); and (f) NCEP surface zonal current (ms" 1 ). The regression maps of 
NSIPP (NCEP) fields are computed using the NSIPP (NCEP) Nino3 index, 
respectively. 
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Figure 9. The leading EOFs of the BY SST and Z20 perturbations derived from the 
NSEPP and NCEP/CFS03 GGCMs. (a) EOF1 of NSIPP BV SST; (b) EOF1 of 
NSIPP BV Z20; (c) EOF2 of NSIPP BV Z20; (d) EOF1 of NCEP BV SST; (e) 
EOF1 of NCEP BV Z20; and (f) EOF2 of NCEP BV Z20. The scale is arbitrary. 














